set.seed(5168)
## responsibility attribution ##

	library(foreign)
	library(mnormt)
	library(ordinal)
  library(MASS)

  data <- read.dta('Netherlands (2012).dta')

  data <- subset(data, data$Party != 1) 
  data <- subset(data, data$Party != 2) 
  data <- subset(data, data$Party != 4) 
  data <- subset(data, data$Party != 5) 
  data <- subset(data, data$Party != 6) 
  data <- subset(data, data$Party != 7) 
  data <- subset(data, data$Party != 8) 
  data <- subset(data, data$Party != 9) 
  data <- subset(data, data$Party != 10) 

  model <- polr(as.factor(Resp_attribution) ~ Prime_minister + Cabinet + Opposition + Perc_seats + Prime_minister * Perc_seats + Cabinet * Perc_seats + Opposition*Perc_seats, data = data, Hess=TRUE, method = 'probit')
	summary(model)

	pars <- as.matrix(model$coefficients)
  pars[8:11] <- as.matrix(model$zeta)     #since this is using a different package I need to add the intercepts to hte vector that takes all of the coefficients
	vce <- vcov(model)


	simCoef <- rmnorm(1000, as.vector(pars),as.matrix(vce))

	summary(model)
	summary(simCoef)

	pmHi <- c()
	pmLo <- c()
	pmMed <- c()

  partnerHi <- c()
  partnerLo <- c()
  partnerMed <- c()

	supHi <- c()
	supLo <- c()
	supMed <- c()

	oppHi <- c()
	oppLo <- c()
	oppMed <- c()

  otherHi <- c()
  otherLo <- c()
  otherMed <- c()

  part_pred <- c()
  opp_pred <- c()
  
  marg_effect <- c()
  margeff_Hi <- c()
  margeff_Lo <- c()
  margeff_Med <- c()

	for(i in 0:25){
	## partner
	agg <- simCoef[, 1] * 0 + 
	  simCoef[, 2] * 1 + 
	  simCoef[, 3] * 0 + 
	  simCoef[, 4] * i + 
	  simCoef[, 5] * 0 + 
	  simCoef[, 6] * i + 
	  simCoef[, 7] * 0  
		
		pred1 <- pnorm(simCoef[, 1], agg, sd = 1)
		pred2 <- pnorm(simCoef[, 2], agg, sd = 1) - pred1
		pred3 <- pnorm(simCoef[, 3], agg, sd = 1) - pred2 - pred1
		pred4 <- pnorm(simCoef[, 4], agg, sd = 1) - pred3 - pred2 - pred1
		pred5 <- 1 - pred1 - pred2 - pred3 - pred4
		
		part_pred <- pred5 * 5 + 
			pred4 * 4 + 
			pred3 * 3 + 
			pred2 * 2 + 
			pred1

		partnerHi[i + 1] <- quantile(pred, 0.975)
		partnerLo[i + 1] <- quantile(pred, 0.025)
		partnerMed[i + 1] <- quantile(pred, 0.5)

	## opposition
	agg <- simCoef[, 1] * 0 + 
	  simCoef[, 2] * 0 + 
	  simCoef[, 3] * 1 + 
	  simCoef[, 4] * i + 
	  simCoef[, 5] * 0 + 
	  simCoef[, 6] * 0 + 
	  simCoef[, 7] * i 
		
		pred1 <- pnorm(simCoef[, 1], agg, sd = 1)
		pred2 <- pnorm(simCoef[, 2], agg, sd = 1) - pred1
		pred3 <- pnorm(simCoef[, 3], agg, sd = 1) - pred2 - pred1
		pred4 <- pnorm(simCoef[, 4], agg, sd = 1) - pred3 - pred2 - pred1
		pred5 <- 1 - pred1 - pred2 - pred3 - pred4
		
	 opp_pred <- pred5 * 5 + 
			pred4 * 4 + 
			pred3 * 3 + 
			pred2 * 2 + 
			pred1

		oppHi[i + 1] <- quantile(pred, 0.975)
	  oppLo[i + 1] <- quantile(pred, 0.025)
	  oppMed[i + 1] <- quantile(pred, 0.5)

  	marg_effect = part_pred - opp_pred
  	margeff_Hi[i + 1] = quantile(marg_effect, 0.975)
  	margeff_Lo[i + 1] = quantile(marg_effect, 0.025)
  	margeff_Med[i + 1] = quantile(marg_effect, 0.5)
    }

#seats sequence
seats <- seq(0, 25, 1)

#Now create the figure
par(mar=c(3.5,6.2,1.5,1.5))    
plot(1, 1, type = 'n',
     xlim = c(0, 25),
     ylim = c(0, 1.5),
     xlab = '',
     ylab = '',
     main = '', 
     axes = FALSE)
axis(1,at=seq(0,25,5),labels=T)         
axis(2,at=seq(0,1.5,0.50),labels=T)
polygon(c(seats, rev(seats)), c(margeff_Hi, rev(margeff_Lo)), col = rgb(red = 0.1, 0.1, 0, 0.2), border = FALSE)
lines(seats, margeff_Med, col = rgb(0.1, 0.1, 0, 1), lwd = 3)
abline(h=0,col="red")

mtext("Change in \nresponsibility attribution", side=2, line=2.5, cex=1.2)       
mtext("Perceived seat %", side=1, line=2.2, cex=1.2)                 

dev.copy(png,'Figure 8b NL 2012 - Marginal effect of support party categorization.png', height = 490, width = 633)  
dev.off()  
